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Supplementary Information 
1. Computational Methods 

1.1. Automated extraction and reduction of prior information from 
signaling databases 

Requirements for prior information. There is a growing number of signaling 
databases (Bader et al., 2006) that capture pathway information in high detail. 
Compared to the interaction networks with relatively low level of detail, pathway 
databases contain information about the phosphorylation states of proteins and 
interactions specific to these states. Although this corpus is highly valuable for 
phospho-proteomic analysis, it remains mostly untapped. To utilize this 
knowledge, one needs algorithms to map the experimentally measured proteins 
and phospho-proteins to their pathway database equivalents, find connections 
between them and reduce the outcome to a format that can be used for inference 
and analysis. 

The Algorithm We developed an algorithm and a software tool (Prior Extraction 
and Reduction Algorithm, PERA) to automatically extract prior information from 
multiple signaling databases in the BioPAX (Demir et al., 2010) format and 
generate a prior information network. PERA takes a list of (phospho)proteins 
identified by their HGNC symbols (e.g. AKT1), phosphorylation sites (e.g. pS473) 
and their molecular status {i.e., activating or inhibitory phosphorylation, total 
concentration) as input and then finds directed signaling paths between these 
entities. These paths are then reduced to directed interactions between signaling 
molecules represented in a Simple Interaction Format (SIF). 

The PERA framework. The prior information is extracted from components of 
the Pathway Commons 2 database in three steps: First, using the paths-between 
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graph query algorithm (Dogrusoz et al., 2009), PERA generates a sub-graph 
(i.e., graph-of-interest) of Pathway Commons, which contains all the input 
proteins and all known connections within their first neighborhoods. Second, 
using the phosphorylation and activity state information, input entities are 
mapped to the corresponding protein states in the graph-of-interest. During this 
mapping step, protein states that do not match with either the corresponding 
annotation for phosphorylation or activity state are filtered out. Phosphorylation 
site mismatches up to 6 residues are tolerated during the filtering step to account 
for phosphorylation site ambiguities due to either database curation errors or 
cross-organism annotations. Third, paths that result in the addition or removal of 
a phopho-group at a phosphorylation site are extracted and mapped to 
phosphoprotein nodes. For total protein nodes, all non-phosphoprotein specific, 
directed signaling paths are included. For this application, the maximum 
allowable graph-query distance limit was set to 1 and only the Reactome 
(Matthews et al., 2009) and NCI-Nature PID pathway (Schaefer et al., 2009) data 
resources were used. Although we limited ourselves to short path distances and 
two pathway databases that we were most familiar with, PERA can be applied to 
extract information from any pathway database that exports to BioPAX and can 
be configured for searching paths of arbitrary length. The PERA software is 
available at http://bit.ly/bp_prior as free software under LGPL 3.0. 

Key considerations for extracting pathway information: Two major semantic 
issues need to be taken into account when extracting state-specific pathway 
information from databases. First semantic issue stems from the ambiguities in 
mapping proteins with multiple phosphorylation sites. When a protein with 
multiple phosphorylation sites is experimentally profiled with an antibody, which 
recognizes only a single phosphorylation site, the antibody will actually bind to a 
heterogeneous mixture of phospho-states provided that the epitope is 
phosphorylated (e.g., anti-AKTpS473 Ab may bind to both AKTpS473 and 
AKTpS473_pT308 but not to AKTpT308). For proteins with multiple observed 
phosphorylation sites, this might lead to semantic conflicts since a double 
phosphorylated node should be mapped to both observations (i.e. single and 
double phosphorylated states). We included an optional "strict" mapping scheme 
to map only the phosphoproteins that exactly match the observed node - in our 
case always single phosphoproteins (e.g., the epitope of anti-AKTpS473 Ab is 
mapped to AKTpS473 but not to AKTpS473_pT308). Since our inference 
algorithm is much more tolerant of missing interactions compared to false 
interactions, we opted to use this flag for this application. 

The second semantic issue stems from the fact that pathway databases are often 
curated from multiple independent studies, spanning multiple cellular states, cell 
and tissue types and even multiple model organisms. As a result, they are a 
superimposition of possible interactions over a wide range of spatiotemporal and 
genetic contexts. On the other hand databases cover only a subset of all possible 
contexts. In our case, we expect only a subset of the interactions in the pathway 
databases to be active in our cell lines and cover only a subset of the interactions 
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that we observe. This observation necessitates incorporating prior information as 
"soft" restraints for network inference. For this purpose, we devised a modified 
cost function, which includes a term for prior information (See below. Also see 
(Molinelli, Korkut, Wang et al., 2013)). 

1.2. Modified cost function for network inference with prior 
information 

In order to generate predictive network models of signaling, we have quantified 
the cost of W by an objective cost function C(W) that penalizes (a) discrepancies 
between predicted ( xf ) and experimental ( xf* ) values of the system variables 
at a time points {ti} in condition u, and (b) the number of non-zero interactions 
with an L0 norm (Molinelli, Korkut, Wang et al., 2013; Nelander et al., 2008). 
Here we have also incorporated prior information to construct network models 
with higher predictive power even if the models are constructed with limited 
experimental data. The newly introduced third term in the cost function accounts 
for the prize introduced when the inferred wy conforms to the prior information 
(Wij Pnor ). The modified cost function with prior information term is formulated as in 
equation S1. 

Equation S1. Modified Cost function 

C(W) = C SSE (W) + c complexity (W) + C prior (W) 

LNM NN NN 

c(w) = (3£X£( x i( t 1 )- x i*( t .)) +^IIS(w ll )+IIn(w, ; ) 

1 i 11 i }*i i j*i 
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In Equation S1, the first term accounts for the discrepancies between predicted ( 
xf ) and experimental ( xf* ) values of the system variables at a time point {ti} in 
condition u, for a particular network configuration ( W= {wy} ). The second term is 
the complexity factor with an L0 norm, which reduces the number of nonzero 
interactions in a network configuration and ensures that resulting network models 
are sparse (Molinelli et al., 2013). The final term, n(W) is the cumulative prior 
cost function for W={wy}. n( w y = w) has a negative real value and reduces the 
overall model cost if the inferred u> conforms to the prior information. Here we 
will formulate the newly introduced prior information term in the modified cost 
function. 

Generalized prior information prize function. In order to derive the prior 
information prize (PIP) function, we first assume the value of inferred wy is 
normally distributed around an expected prior value, E[wy pnor ] with a standard 
deviation, oy pnor . The prior probability model for each interaction in the prior 
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network is formulized using the normal distribution probability density function 
around E[Wij pnor ] as in equation S2. 

Equation S2. Probability model of prior observations 

P pnor (w..) = . exp(-^ . \ ) 

IJ af'V^ 2(a P i nor ) 2 

The equation S2 implies that the probability of inferring a particular wy value and 
consequently, the prize introduced to the cost function increase as wy 
approaches E[wy pnor ]. We calculate the prize for the fit between the inferred wy 
and the E[Wy prior ] with an inverse Boltzmann conversion of the probability model 
(Jaynes, 1957). 



Equation S3. Error model of individual prior observations 
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The first term introduces the maximum possible prize for each interaction that is 
represented in the prior information network. Second term penalizes the 
discrepancies between the wy and the E[wy pnor ]. The penalty for the discrepancy 
is zero when wy = E[Wij Pnor ] and the prior prize, n( w y) assumes the highest 
possible value. Ky(ln(1/ (Oij Pnor V27r)) is a constant analogous to the inverse 
temperature in statistical physics. The constant Ky and the resulting n(wy) is zero 
if the interaction between nodes i and j is not contained in the prior network 
model. By summing over all i and j, the cumulative prior prize for all W={wy} 
becomes: 



Equation S4. Cumulative error model of prior observations 
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Even though the current PIP function is derived with the normal distribution 
assumption, a variety of alternative distributions such as bimodal distribution can 
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be fitted in Equation S2. Such efforts will lead to a variety of custom designed 
PIP functions when necessary. 

Estimation of E[Wij prior ] and Oij prior 

An important requirement for the implementation of equation S4 is the accurate 
estimation of E[Wij pnor ] and oy prior . Here, we propose two estimation strategies that 
are alternative to systematic extraction of priors from databases. 

Priors from experiments. Our first strategy relies on use of biochemical 
measurements with multiple biological replicates to quantify the influence of node 
j on node I (wy). The E(wy pnor ) and oy pnor can be estimated from the mean value 
(<wy pnor >) and standard deviation (oy) of the readouts from biological replicates 
in such measurements. Depending on the nature of the interaction under study, a 
variety of different biochemical methods can be utilized. These methods include 
but are not limited to enzyme activity assays, quantification of protein-protein 
interactions (Selbach and Mann, 2006), or additional RPPA readouts under 
specific perturbation conditions. 

Priors from network models. As we build network models of signaling in 
diverse biological contexts using perturbation data and inference algorithms, a 
vast amount of information on signaling properties will become available. 
Consequently, we will be able to estimate the E(wy pnor ) and Oy pnor for newly 
inferred models from available network models. Previously inferred {wy} values 
and standard deviation in BP-generated probability distributions will serve as a 
basis to estimate the E(wy pnor ) and oy pnor , respectively. 

The simplified PIP function 

The current form of the prior information is limited to a set of binary interactions 
due to the qualitative nature of the databases, from which we extract the 
information. For instance, the prior information stored in the databases may 
correspond to activating (wy>0), inhibitory (wy<0), or generic (wy^O) interactions. 
In such situations, wy prior has a uniform distribution within the defined boundary. 
Thus, the E[wy prior ] can be set as equal to any value within the interval, [wy min " prior , 
Wj .max-pnoi-j p Qr exam p| e a p r j or f or a positive interaction between nodes i and j is 

represented by an interval of (0, wy max ], where wy max is the maximum allowed 
edge strength in the inference scheme. Consequently, the squared distance error 
in equation S3 vanishes and the prior prize assumes a fixed value over the 
defined prior interval. 
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Equation S5. Simplified error model of prior information 
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The binary nature of the interactions also implies a generic k value 
(K=Kij(ln(1/(Oij pnor V27r)), icy >0, V i,j £ N ) for all interactions represented in prior 
information network and k=0 for all other situations. The resulting prior 
information term in the cost function is a step function as shown in the methods 
function. 

Equation S6. Simplified Error model of individual prior observations 

TlCWj- # 0) = -K and T|(w y = 0) = 0 Generic prior information (w^O) 

TltWy > 0) = -K and T^w- < 0) = 0 Prior information for an activating interaction (w t . > 0) 
tl(Wij < 0) = -K and ri(Wij > 0) - 0 Prior information for an inhibitory interaction (w» < 0) 
TiCWy) = 0 No Prior exists for w ;j 

For each interaction in the prior information network, n( w y) has a negative value. 
Thus, the prior information introduces a prize in the overall cost function. The 
formulation and implementation into the BP equations (see main text) create a 
soft-restraint for prior information extracted from signaling databases. In our 
probabilistic inference framework, a prior information is accepted as an edge 
(wy^O) only if the prize introduced can override any discrepancy between 
predicted (Xi M ) and experimental ( Xi M * ) values of the system variables. Thus the 
context specific character of inferred models is preserved when the cost function 
in equation S1 is used. 

1.3. Network inference 

The computational/experimental network modeling procedure is described in 
Figure 1. The theoretical details of the BP-based network inference is described 
elsewhere (Molinelli, Korkut, Wang et al. 2013). The BP-guided decimation 
algorithm is described in Figure S1. The mean-field approximation in belief 
propagation inference algorithm is summarized below. 

Gaussian integration of cavity parameter update. Information from a 
sufficiently large set of non-cavity constraints and parameters are integrated to 
update the cavity parameters. As a corollary to the central limit theorem, we can 
assume that the mean-field effect of non-cavity parameters on the cavity update 
follows a Gaussian distribution. The mean and variance of the Gaussian 
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distribution are the means and variances of the global distributions given by 
individual distributions of non-cavity P^(wy) vj^k. Thus, we replace the sum over 
multivariate configurations (Equation 2.B in main text) of all non-cavity 
parameters with a single Gaussian integration over the mean-field of non-cavity 
configurations. 

Equation S7. Gaussian approximation to the cavity parameter update 
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Note that, BP parameters (temperature factor, p=2, complexity factor, A=5) are 
adapted for all computations after a systematical analysis of both parameters 
described in Molinelli, Korkut, Wang et al. Prior information parameter, k is set as 
(k= A =5) to ensure prior information and sparsity restraints in the models balance 
each other. When prior information is used to construct the models, the 
complexity term, A is set to 2.5 for phenotypic nodes in order to approximately 
match the complexity around the phenotypic nodes and rest of the networks. 
When the parameter space is sufficiently large, this approximation significantly 
reduces the computational complexity of the problem without loss of overall 
accuracy. By means of the iterative approach with the mean-field Gaussian 
approximation, the time-complexity of the problem is strongly reduced and the 
obstacles imposed by combinatorial complexity are circumvented. 



2. Experimental methods 

Cell cultures and Perturbation experiments. In perturbation experiments, the 
drugs are applied to cell cultures after SkMel133 cells are grown to about 40% 
confluence in complete RPMI-1640 medium (10% heat-inactivated fetal bovine 
serum, 100 units/ml each of penicillin and streptomycin, and incubated at 37°C in 
5% CO2) in 6-well plates. 24 hours after drug administration, the perturbed cells 
are harvested. In control experiments (i.e. no-drug condition), cells are treated 
with the DMSO drug vehicle for 24 hours. 

Reverse Phase Protein Arrays. Proteomic response profiles are measured 
using reverse phase protein arrays (MD Anderson Cancer Center RPPA Core 
Facility) (Tibes et al., 2006). For sample preparation, cell pellets are lysed in 
RPPA buffer (1% Triton X-100, 50mM HEPES pH7.4, 150mM NaCI, 1.5mM 
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MgCI2 , 1mM EGTA, 100mM NaF, 10mM Na pyrophosphate, 1mM Na3V04, 
10% glycerol, freshly added protease and phosphatase inhibitors). After cells are 
lysed, total protein concentrations in cell lysates are measured with the Bradford 
assay and the final protein concentrations are adjusted to 1-1 .5mg/ml. Samples 
are boiled in SDS sample buffer without bromophenol blue for 5 minutes at 95°C. 
Cell lysates are spotted on nitrocellulose-coated slides as described in (Tibes et 
al., 2006). Antibody staining intensities are quantified using the MicroVigene 
automated RPPA module (VigeneTech, Inc.) and the standard RPPA protein 
concentration normalization procedure (Neeley et al., 2009) is followed. The 
proteomic readouts are log normalized with respect to the corresponding 
untreated condition readouts. We have eliminated those readouts with intra or 
inter slide coefficient of variation (CV) > 0.15 (I.e., low reproducibility) and low 
degree of staining by antibodies. 100 proteomic entities are chosen for further 
analysis. Next, we excluded the proteomic entities that do not respond to at least 
a single perturbation condition from the network models. This has lead to highly 
reliable response measurements on a final set of 82 proteomic entities, which 
have been used in the network models. 

Quantitative phenotypic assays. All phenotypic measurements are made in 
perturbation conditions identical to those in proteomic measurements. Cell 
viability is measured 72 hours after perturbation using the Resazurin assay. Each 
well is treated with the Resazurin (Sigma-Aldrich, Catalog # R7017) at a final 
concentration of 44 uM for 1 hour. The fluorescent signals are measured at 530 
nm excitation wavelength and 590 nm emission wavelength. Standard curves of 
cell numbers are used to calculate the cell numbers in each sample. The cell 
cycle progression phenotypes are measured with flow cytometry analysis. Cells 
are seeded in 6-well plates. 72 hours after drug application, cells are harvested 
by trypsinization, including both suspended and adherent fractions, and washed 
in cold PBS. Cell nuclei are prepared as described in (Nusse and Kramer, 1984) 
and cell cycle distribution was determined by flow cytometric analysis of DNA 
content using red fluorescence of 488 nm excited ethidium bromide-stained 
nuclei. The percentage of cells in the G1, G2/M, and S phases and sub-G1 
fraction are recorded. 
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Supplementary Figures: 

BP-guided decimation algorithm 
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Figure S1 BP-guided decimation algorithm is used to construct executable, 
individual network model solutions from BP generated probability distributions for 
each edge strength value (wij). The algorithm chart depicts one round of BP- 
guided decimation to generate a single model solution. Consecutive runs of BP- 
guided decimation algorithm leads to construction of a network model solution 
ensemble. 
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Continued 



Proteomic response to targeted drugs in SkMel133 cells -continued 
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Figure S2 Clustering and pathway analysis of proteomics data. We 

characterized the proteomic signatures in the response data with a pathway 
analysis guided by hierarchical clustering. A. The two-way clustering analysis of 
the experimental response map reveals distinct proteomic signatures of response 
to drugs targeting different signaling pathways (See Figure 2 in main text). The 
Cluster software is used for the two way hierarchical clustering with correlation- 
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based distance metric and average-linkage method. We identified 3 major (C1, 
C2 and C4) and 3 minor (C3, C5, C6) proteomic response clusters. B. We 
extracted the signaling interactions between the proteomic entities that fell into 
each cluster. The signaling interactions within each proteomic cluster are 
extracted from Pathway Commons 2 database using the PERA algorithm. We 
simplified the pathway diagrams by removing the post-translational modifications 
and merging the nodes associated with identical genes (e.g., AKTpT308, 
AKTpS474 and AKT are merged to AKT node). The resulting diagrams displayed 
the known pathways associated with proteomic clusters, whose members gave 
similar response to targeted agents. The cluster guided pathway analysis 
suggests that functionally related proteins (i.e. proteins on the same or related 
pathways) are enriched in distinct clusters. The phospho-proteomic entities on 
MAPK and PI3K/AKT pathways are enriched in cluster 1 (C1) of the response 
map. The total protein measurements related to MAPK and PI3K/AKT pathways 
are enriched in a related but distinct cluster (C2). The level of various apoptotic 
proteins increase in response to most targeted drugs and form another distinct 
cluster (C4). Interestingly, EGFR (Both phospho and total level) and a set of 
EGFR related docking proteins (e.g., SHC, 14-3-3-P) are also enriched in C4 
possibly due to the increase in the expression of those entities in response to 
targeted drugs such as MEKi and PI3KL The observed increase in EGFR level 
suggests a potential feedback loop that emerges from downstream elements of 
MAPK and PI3K/AKT pathways, which are targeted by multiple agents in this 
study. C. Proteomic response to single agent perturbations. Proteomic levels in 
SkMel133 cells change in response to targeted drugs as measured with RPPA. 
For each drug, we listed the top 10 decreasing and increasing phospho-proteins 
out of 138 proteomic entities profiled with RPPA (Figure S1A). The protein levels 
are normalized to no drug condition and given in linear space. 
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A Qualitative prior model 
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Statistical properties of the solution ensemble 
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Figure S3. A. The prior model of signaling. The prior information model is 
extracted from Reactome and NCI PID using the PERA algorithm (See methods 
and supplementary methods for PERA algorithm and use of the prior model in 
inference). The prior model served as a bias in BP-based network inference. 
Green arrows represent priors for activating edges, red arrows represent 
inhibitory edges and black arrows represent generic edges (i.e. activating or 
inhibitory). B. Distribution of edges in the solution ensemble. In order to 
model cellular signaling and predict response, we computed 4000 signaling 
models and generated a solution ensemble (See main text). A. Edge frequencies 
(f(wy)) in the solution ensemble. The y-axis represents the frequency values for 
nonzero edge strengths (f(|wy|>0.2)) in the ensemble. The X-axis represents the 
interactions ordered by their f(|wy|>0.2). The long tail that corresponds to the 
edges with frequencies less than 0.05 is not displayed. A set of well reported 
interactions such as the coupling of AKT activity node to phospho-AKT and 
phospho-MEK to phosho-MAPK (ERK) are inferred with |wy|> 0.2 in >99% of the 
model solutions (left). Few examples of the rarely captured signaling interactions 
are given on right. B. The frequency distribution of nonzero edges (|wy|>0.2) has 
a bimodal character. A set of interactions with nonzero edge values is shared by 
nearly all of the inferred network models (0.8 < f( |wy|> 0.2 ) < 1.0). Such 
frequent edges form a stable network skeleton shared by majority of solutions. 
On the other hand, a set of possible interactions have non-zero edge strength 
(wy) values in around 50% of the solutions. Such interactions vary among 
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different model solutions and possibly account for the variations in system 
response predicted with different model solutions. The interactions in the prior 
model are particularly enriched in these two groups. A large fraction of potential 
interactions have very low (0.05< f(|wy|>0.2)<0.20) or near zero (f(|wy|>0.2) < 
0.05) frequencies in the model ensemble. C. Comparison of random vs. actual 
prior information. We tested whether the database driven prior model 
conformed to the experimental data and prior information does not overly restrain 
the inferred models. For this purpose, we constructed and compared network 
models using the pathway driven and randomly generated prior restraints in BP- 
based modeling. We have tested 500 unique, randomly generated prior restraints 
each containing 154 interactions, equal to the number of interactions in the 
pathway driven prior model. For model construction, we first computed the 
probability distribution of edge values (P(Wy)) for each possible interaction. Next, 
we constructed the models by assigning the edge value (Wy) for which BP- 
generated P(Wy) is maximum. We compared the models generated with different 
prior models for their prior scores (i.e. number of edges, which were accepted by 
the algorithm and also contained in the prior model). The models generated 
using the database driven priors had significantly higher prior scores (ps) 
compared to randomly generated models (p < 0.05 Student's t-test for HO: 

Ups rand ° m = ^database-drive^ NQte ^ W0 performed 500 parallel BP TUnS With 

identical database driven prior models and experimental data. We observed a 
low degree of variation in the prior scores for the models with database driven 
priors. This is indicative of some solution degeneracy in the models inferred with 
database driven priors. In predictive network modeling, we accounted for the 
solution degeneracy by re-computing the P(Wy) at the beginning of each BP- 
guided decimation step (See main text and methods) and ensuring that the final 
solution ensemble contains models that sample every possible degenerate 
solution. 
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Cell viability 
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Figure S4 A. Prediction of phenotypic responses to in silico, combinatorial 
perturbations. We computationally predicted cell viability and cell cycle arrest 
(G1, S, G2 and G2M) in response to combinatorial perturbations of all proteomic 
nodes in network models (See Table 1 and Figure 5-6 in main text). In silico 
perturbations include paired combinations of 94 perturbations in 4 different 
concentrations (ICO-20-40-60-80 with respect to the basal level of a particular 
node) leading to more than 70000 unique perturbation conditions. Each in silico 
perturbation is applied to 4000 inferred network models. This leads to a 
simulation of -28 million unique conditions for 99 different readouts, a throughput 
level inaccessible by experimental screening methods. The heatmaps display the 
phenotypic response landscape to combinatorial perturbations. Each data point 
on heatmaps is the response to a perturbation averaged over predictions from 
4000 network models A. Cell viability response to perturbations. The desired 
response from a perturbation is reduction in cell viability. (Red: decreased cell 
viability, green: increased cell viability). (B-E). Same as in A but for cell cycle 
progression phenotypes. For cell cycle phenotypes, the desired response is an 
increase in the cell cycle arrest (Red: increased cell cycle arrest. Green: 
Decreased cell cycle arrest). F. The top ten most effective single-agent in silico 
perturbations for each phenotype. The listed proteomic entities for each 
phenotype are potential drug targets for slowing down the growth of melanoma 
cell. The error bars represent the ±SEM over the simulations of all model 
solutions. For complete prediction heatmaps for phenotypes see Supplementary 
figure 6-7. 
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D- Uncropped gels 
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Figure S5. A. The cell viability response to CDK4i and PKCi. The IC50 for 
CDK4i is ~2[iM in SkMel133 cells. The maximum response is reached ~10^iM B. 
The IC50 for CDK4i is -l^iM and maximum response is reached at 4[iM in 
SkMel133. CDK4i (2\i) and PKCi (3.5-7|xM) were used in original input 
perturbation data. The models recapitulated the effect of high dose PKCi and 
CDK4i to inhibit cell viability. However, no significant effect on viability is 
observed in the nM range. B. Changes in c-Myc level in response to 
perturbations in Skmel133 cells as measured in RPPA experiments. Each 
data point is log normalized with respect to the c-Myc level in unperturbed 
condition. c-Myc level is highest in the presence of CDK4i. Various drug 
combinations that include STAT3 and mTOR inhibitors follow CDK4i 
combinations. cMyc level is lowest when cells are perturbed with MEK, BRAF, 
SRC, PKC and HDAC inhibitors, respectively. Each data-point is the average of 
RPPA readouts from three replicates. C. Cell cycle arrest response of 
melanoma cells to JQ1 and MEKi combination. 51% and 88% of cells are in 
G1 stage 24 hours after JQ1 and MEKi treatment, respectively. JQ1+MEKi 
combination has an increased effect on G1 cell cycle arrest (G1=92%). 39% of 
cells are in G1 when they are not treated with drugs. Error bars in right panel: 
±SE in 3 biological replicates. D. Uncropped versions of the gels presented in 
Figure 6B. In the fourth gel, the observed band ~100kDa is potentially the BRD4 
B or C isoform. 
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Supplementary Tables: 



Table S1. Drugs used in perturbation experiments 
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0.6 


ZSTK474 


P13K 


AKTpS473 






PD0325901 


MEK1/2 


MAPKlpT202 


1.5 


3 


Stattic 


STAT3 


STAT3pY705 


600 


1200 


HNHA 


HDAC 


P27/Kipl 


6,000 


12,000 




BRAF V600E 








RO-3 1-7549 


PKC 


S6pS240 


3500 


7000 




P6 


JAK 


GSK3pS21 


20 


40 








2,400 


4,800 


Nutlin 


MDM2 


4EBPlpS65 


3,000 


6,000 



25 



Table S2 Proteins that respond to at least one perturbation and exist in models 
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Table S3. Perturbation conditions 
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Dose2 
(nM) 


Index 


Perturbationl 


Dosel 
(nM) 


Perturbation2 


Dose2 
(nM) 


-i 


MFKi 


i c 






At* 








innn 












2 


MEKi 


1.5 


HDACi 


6,000 


47 


BRAFi 


60 


JAKi 


20 






















4 


MEKi 


1.5 


JAKi 


20 


49 


BRAFi 


60 


SRCi 


2400 












50 












6 


MEKi 


1.5 


PKCi 


3500 


51 


BRAFi 


60 


raTORi 


0.3 


7 


MEKi 


1.5 


SRCi 


2400 


52 


PKCi 


3500 






8 


MEKi 


1.5 


STAT3i 


600 


53 


PKCi 


3500 


mTORi 


0.3 



10 

11 

12 



MEKi 
MEKi 
AKTi 



1.5 PI3Ki 



600 



10,000 



55 
56 
57 



CDK4i 
CDK4i 
CDK4i 



2,000 
2,000 
2,000 



14 


AKTi 


5,000 


MEKi 


1.5 


59 


CDK4i 


2,000 


MDM2i 


3000 


15 


AKTi 


5,000 


HDACi 


6000 


60 


CDK4i 


2,000 


JAKi 




16 


AKTi 


5,000 


MDM2i 


3000 


61 


CDK4i 


2,000 


BRAFi 


60 






5,000 


JAKi 


20 


62 


CDK4i 


2,000 


PKCi 




18 


AKTi 


5,000 


BRAFi 


60 


63 


CDK4i 


2,000 


SRCi 


2400 


19 


AKTi 


5,000 


PKCi 


3500 


64 


CDK4i 


2,000 


STAT3i 




20 


AKTi 


5,000 


SRCi 


2400 


65 


CDK4i 


2,000 


mTORi 


0.3 






5,000 


STAT3i 


600 


66 


CDK4i 


2,000 


PI3Ki 




22 


AKTi 


5,000 


mTORi 


0.3 


67 


SRCi 


2,400 






23 


AKTi 


5,000 






68 


SRCi 


2,400 


PKCi 




24 


HDACi 


12,000 






69 


SRCi 


2,400 


mTORi 


0.3 




HDACi 


6,000 






70 


SRCi 


4,800 






26 


HDACi 


6,000 


MDM2i 


3000 


71 


STAT3i 


600 






" 


HDACi 


6,000 


PKCi 


3500 


72 


STAT3i 


600 


HDACi 




28 


HDACi 


6,000 


SRCi 


2400 


73 


STAT3i 


600 


MDM2i 


3000 




HDACi 


6,000 


mTORi 


0.3 


74 


STAT3i 


600 


PKCi 




30 


MDM2i 


3,000 






75 


STAT3i 


600 


SRCi 


2400 




76 STAT3i 600 mTORi 0.3 


32 


MDM2i 


3000 


SRCi 


2400 


77 


STAT3i 


1200 






33 MDM2i 3000 mTORi 0.3 


78 


mTORi 








34 


MDM2i 


6,000 






79 


mTORi 


0.6 










36 


JAKi 


20 


HDACi 


6000 


81 


PI3Ki 


600 


HDACi 


6000 




JAKi 


20 


MDM2i 


3000 


82 


PI3Ki 


600 


MDM2i 




38 


JAKi 


20 


PKCi 


3500 


83 


PI3Ki 


600 


JAKi 


20 






40 


JAKi 


20 


STAT3i 


600 


85 


PI3Ki 


600 


PKCi 


3500 




42 


JAKi 


40 






87 


PI3Ki 


600 


STAT3i 


600 






44 


BRAFi 


60 






89 


PI3Ki 


1,200 
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Table S4. Predicted G1 -arrest response to combined targeting of c-Myc and 
other nodes (Perturbations on c-Myc + X). The response values are in log2 and 
normalized with respect to the response to single agent c-Myc perturbation (G1 C - 
Myc+x/G1 c-Myc)- The top 20 paired perturbation conditions are ranked based on the 
response to the lowest dose (c-Myc_IC20). The _ICX tag determines the 
strength of the in silico perturbation. The X (e.g, 20, 40) represents the % 
decrease in the target node activity at the given perturbation strength. 





c-Myc_IC20 c- 


Myc_IC40 


c-Myc_IC60 


c-Myc_IC80 


aMEK_IC80 




1.971 


1.870 


1.835 


aMEK_IC60 


2.253 


1.924 


1.824 


1.791 


aMEK_IC40 


2.071 


1.786 


1.698 


1.668 


CyclinDl_IC80 


1.928 


1.699 


1.630 


1.607 


CyclinDl_IC60 1.893 1.670 1.606 1.581 


aBRAFm_IC80 


1.889 


1.677 


1.611 


1.586 


aBRAFm_IC60 


1.845 


1.640 


1.575 


1.554 


CyclinDl_IC40 


1.795 


1.594 


1.536 


1.515 


MAPKpT202_IC80 


1.725 


1.535 


1.478 


1.458 


aBRAFm_IC40 


1.721 


1.540 


1.482 


1.463 




MAPKpT202_IC60 


1.686 


1.506 


1.453 


1.434 




PLK1IC80 


1.679 


1.531 


1.479 


1.462 


aSRC_IC60 


1.678 


1.526 


1.474 


1.456 


PLK1IC60 


1.670 


1.523 


1.470 


1.452 












PLK1IC40 


1.639 


1.493 


1.445 


1.430 
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